---
title: "Variant Calling"
author: "Pierre Barry"
date: "`r format(Sys.time(), '%d %B, %Y, %H:%M')`"
output: 
  flexdashboard::flex_dashboard:
    theme: paper
    orientation: rows
    social: menu
    source_code: embed
    vertical_layout: scroll
---

```{r global, include=FALSE}
list.of.packages <- c("RColorBrewer","stringr","kableExtra","formattable","ggplot2","plotly","DT")

for (i in list.of.packages){
  if (i %in% installed.packages()[,"Package"] == FALSE){
    install.packages(i)
  }
  eval(bquote(library(.(i))))
}
```

```{r, include = FALSE}
color_med_atl=data.frame(Location=c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                         Col=brewer.pal(n = 4, name = "RdBu"))
```


VCF {data-icon="ion-ios-settings"}
=======================================================================

```{r}
load(file=paste(path,"/output/variant_calling/VCF_SPECIES.Rdata",sep=""))
load(file=paste(path,"/output/variant_calling/VCF_STATS_INDIV.Rdata",sep=""))
load(file=paste(path,"/output/variant_calling/VCF_DEPTH_QUALITY_PER_SITE.Rdata",sep=""))
LOCATION=c()
for (i in 1:nrow(VCF_STATS_INDIV)){
  LOCATION[i]=location=str_sub(VCF_STATS_INDIV$SAMPLES[i],start=6,end=7)
}
VCF_STATS_INDIV$LOCATION=LOCATION

for (j in 1:nrow(VCF_SPECIES[[3]])){
  VCF_SPECIES[[3]][j,]$COUNT=VCF_SPECIES[[3]][j,]$COUNT/
    (VCF_SPECIES[[1]]$NBR_SNPS[which(VCF_SPECIES[[1]]$SPECIES==VCF_SPECIES[[3]]$SPECIES[j])])
}

labels_species=c()

for (sp in 1:length(levels(factor(VCF_STATS_INDIV$SPECIES)))){
  
  labels_species=c(labels_species,
                   paste(levels(factor(VCF_STATS_INDIV$SPECIES))[sp]," \n n = ",table(VCF_STATS_INDIV$SPECIES)[sp],sep=" "))
  
}

labels_species=as.vector(labels_species)

VCF_SPECIES$LOCATION=factor(VCF_SPECIES$LOCATION,levels=c("Li","Mu","Fa","Ga"))
VCF_STATS_INDIV$LOCATION=factor(VCF_STATS_INDIV$LOCATION,levels=c("Li","Mu","Fa","Ga"))
```

Row {.tabset}
-----------------------------------------------------------------------

### Filter 

```{r}
filter_variant = read.csv(paste(path,"/output/variant_calling/variant_filter.csv",sep=""))
filter_variant = filter_variant[,-1]

filter_variant$PERCE_NBR_SNPS_PASS_FROMINDEL5BP <- round(filter_variant$NBR_SNPS_PASS*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$PERCE_NBR_SNPS_PASS_FROMINDEL5BP <- color_bar("lightgreen")(filter_variant$PERCE_NBR_SNPS_PASS_FROMINDEL5BP)

filter_variant$PERCE_NBR_SNPS_PASS_FROMSNPRAW <- round(filter_variant$NBR_SNPS_PASS*100 / filter_variant$NBR_SNPS_RAW,4)
filter_variant$PERCE_NBR_SNPS_PASS_FROMSNPRAW <- color_bar("lightgreen")(filter_variant$PERCE_NBR_SNPS_PASS_FROMSNPRAW)

filter_variant$NBR_SNPS_PASS <- color_bar("lightgreen")(filter_variant$NBR_SNPS_PASS)

filter_variant$NBR_SNPS_FILTER_ReadPosRankSumm8 <- round(filter_variant$NBR_SNPS_FILTER_ReadPosRankSumm8*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_ReadPosRankSumm8 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_ReadPosRankSumm8)

filter_variant$NBR_SNPS_FILTER_MQRanksSumm12p5 <- round(filter_variant$NBR_SNPS_FILTER_MQRanksSumm12p5*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_MQRanksSumm12p5 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_MQRanksSumm12p5)

filter_variant$NBR_SNPS_FILTER_SOR3 <- round(filter_variant$NBR_SNPS_FILTER_SOR3*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_SOR3 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_SOR3)

filter_variant$NBR_SNPS_FILTER_FS60 <- round(filter_variant$NBR_SNPS_FILTER_FS60*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_FS60 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_FS60)

filter_variant$NBR_SNPS_FILTER_MQ40 <- round(filter_variant$NBR_SNPS_FILTER_MQ40*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_MQ40 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_MQ40)

filter_variant$NBR_SNPS_FILTER_AVGDPcutoff <- round(filter_variant$NBR_SNPS_FILTER_AVGDPcutoff*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_AVGDPcutoff <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_AVGDPcutoff)

filter_variant$NBR_SNPS_FILTER_AVGDP10 <- round(filter_variant$NBR_SNPS_FILTER_AVGDP10*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_AVGDP10 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_AVGDP10)

filter_variant$NBR_SNPS_FILTER_AVGGQ20 <- round(filter_variant$NBR_SNPS_FILTER_AVGGQ20*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_AVGGQ20 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_AVGGQ20)

filter_variant$NBR_SNPS_FILTER_QUAL30 <- round(filter_variant$NBR_SNPS_FILTER_QUAL30*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_QUAL30 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_QUAL30)

filter_variant$NBR_SNPS_FILTER_QD10 <- round(filter_variant$NBR_SNPS_FILTER_QD10*100 / filter_variant$NBR_SNPS_AFTER_INDEL5BP,4)
filter_variant$NBR_SNPS_FILTER_QD10 <- color_bar("lightgreen")(filter_variant$NBR_SNPS_FILTER_QD10)

filter_variant$NBR_SNPS_AFTER_INDEL5BP <- round(filter_variant$NBR_SNPS_AFTER_INDEL5BP*100 / filter_variant$NBR_VARIANT_RAW,4)
filter_variant$NBR_SNPS_AFTER_INDEL5BP <- color_bar("lightgreen")(filter_variant$NBR_SNPS_AFTER_INDEL5BP)

filter_variant$NBR_MULTI_RAW <- color_bar("lightgreen")(filter_variant$NBR_MULTI_RAW)
filter_variant$NBR_INDELS_RAW <- color_bar("lightgreen")(filter_variant$NBR_INDELS_RAW)
filter_variant$NBR_SNPS_RAW <- color_bar("lightgreen")(filter_variant$NBR_SNPS_RAW)
filter_variant$NBR_VARIANT_RAW <- color_bar("lightgreen")(filter_variant$NBR_VARIANT_RAW)

for (j in c(2,3,4,5,18)){
  for (i in 1:length(filter_variant[,j])){
  
    filter_variant[i,j] =  paste(strsplit(filter_variant[i,j],">")[[1]][1],
            ">",
            round(as.numeric(strsplit(strsplit(filter_variant[i,j],"</span>")[[1]],">")[[1]][2])/1e6,3),
            " M",
            "</span>",
            sep=""
      )
      
  }
}


kbl(filter_variant,col.names = c("Species",
                           "Number of variants (raw)",
                           "Number of SNPs (raw)",
                           "Number of indels (raw)",
                           "Number of MNPs (raw)",
                           "SNPs (after removing 5bps around indels, %)",
                           "Filter QD<10 (%)",
                           "Filter QUAL<30 (%)",
                           "Filter AVGGQ<20 (%)",
                           "Filter AVGDP<10 (%)",
                           "Upper depth cutoff",
                           "Filter AVGDPcutoff<XX (%)",
                           "Filter MQ<40 (%)",
                           "Filter FS<60 (%)",
                           "Filter SOR<3 (%)",
                           "Filter MQRankSum<-12.5 (%)",
                           "Filter ReadPosRankSum<-8 (%)",
                           "Number SNPs pass",
                           "SNPs pass (% from indel5bp)",
                           "SNPs pass (% from raw)")
    ,escape = F,align=c(rep('c',times=20))) %>%
  kable_classic ("striped","hover",full_width = T,html_font = "Cambria")
```

### Species

```{r}
rownames(VCF_SPECIES[[1]]) = NULL
VCF_SPECIES[[1]]$SPECIES = c("Atherina boyeri",
                             "Alosa fallax",
                             "Coryphoblennius galerita",
                             "Coris julis",
                             "Dicentrarchus labrax",
                             "Diplodus puntazzo",
                             "Gobius niger",
                             "Hippocampus guttulatus",
                             "Lophius budegassa",
                             "Lithognathus mormyrus",
                             "Merluccius merluccius",
                             "Mullus surmuletus",
                             "Pagellus erythrinus",
                             "Serranus cabrilla",
                             "Spondyliosoma cantharus",
                             "Symphodus cinereus",
                             "Sardina pilchardus",
                             "Sarda sarda",
                             "Syngnathus typhle")
VCF_SPECIES[[1]]$SPECIES = cell_spec(VCF_SPECIES[[1]]$SPECIES,italic = T)
VCF_SPECIES[[1]]$NBR_SNPS <- color_bar("lightgreen")(VCF_SPECIES[[1]]$NBR_SNPS)
VCF_SPECIES[[1]]$TS_TV <- color_bar("lightgreen")(VCF_SPECIES[[1]]$TS_TV)

for (i in 1:length(VCF_SPECIES[[1]]$NBR_SNPS)){
  
VCF_SPECIES[[1]]$NBR_SNPS[i] =  paste(strsplit(VCF_SPECIES[[1]]$NBR_SNPS[i],">")[[1]][1],
        ">",
        round(as.numeric(strsplit(strsplit(VCF_SPECIES[[1]]$NBR_SNPS[i],"</span>")[[1]],">")[[1]][2])/1e6,3),
        " M",
        "</span>",
        sep=""
  )
  
}

kbl(VCF_SPECIES[[1]],col.names = c("Species",
                           "Number of SNPs",
                           "Transitions (Ts)",
                           "Transversions (Tv)",
                           "Ts/Tv")
    ,escape = F,align=c(rep('c',times=5))) %>%
  kable_classic ("striped","hover", full_width = T,html_font = "Cambria")

#kable(VCF_SPECIES[[1]], escape = F) %>%
#  kable_paper("hover", full_width = F)
```

### Substitution 

```{r, fig.height=150}
ppp<-ggplot(VCF_SPECIES[[3]],aes(x=TYPE,y=COUNT))+
  facet_grid(SPECIES ~ .)+
  geom_bar(alpha=0.75,stat='identity',aes(fill=TYPE))+
  scale_fill_viridis_d()+
  theme_classic()+
  xlab("Subsitution type")+
  ylab("Count")

ppp<-ggplotly(ppp,height=2500)  %>% partial_bundle() 
ppp
```

### Table {data-height=750}

```{r}
datatable(VCF_STATS_INDIV, 
          class = 'cell-border stripe',
          rownames = FALSE,
          filter = 'top', options = list(
          pageLength = 100,scrollY=500,autowidth=T,scrollX=100))
```

### Depth

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=DEPTH))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Depth: ",round(DEPTH,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
    scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Depth")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Number of SNPS

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=NBR_SNPS))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Nbre NPS: ",round(NBR_SNPS,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Number of SNPS")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### SNPS same as reference

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=PERCENTAGE_SNPS_SAME_AS_REF))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Percentage same as reference: ",round(PERCENTAGE_SNPS_SAME_AS_REF,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
    scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Percentage same as reference (%)")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Missing genotype

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=PERCENTAGE_MISSING_GENOTYPE))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Percentage of missing genotype: ",round(PERCENTAGE_MISSING_GENOTYPE,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Percentage of missing genotype (%)")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Phased genotype percentage

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=PHASED_GENOTYPE_PERCENTAGE))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Phased genotype percentage: ",round(PHASED_GENOTYPE_PERCENTAGE,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Phased genotype percentage (%)")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Phased genotype number

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=PHASED_GENOTYPE_NBR))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Phased genotype number: ",round(PHASED_GENOTYPE_NBR,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Phased genotype number")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Number of genotype

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=GENOTYPE_NBR))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Nbr of genotypes: ",round(GENOTYPE_NBR,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Nbr of genotypes")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### SNPs Ts/Tv

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=SNP_TS_TV))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Snps Ts/Tv: ",round(SNP_TS_TV,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Snps Ts/Tv")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### SNPs Het/Hom

```{r}
ppp <- ggplot(VCF_STATS_INDIV, aes(x=SPECIES, y=HET_HOM))  +
  geom_violin(fill="chartreuse3",col="black",alpha=0.25) +
  geom_point(aes(col=LOCATION,
                 text=paste("Sample: ",SAMPLES, " \n Snps Het/Hom: ",round(HET_HOM,3),sep="")),size = 1,
             shape=19,
             position='jitter',
             alpha=0.75) +
  stat_summary(fun=median, geom="point", size=2, color="orange")+
  scale_x_discrete(name="Species",
                   labels=labels_species)+
  scale_colour_manual(name = "Location",
                    labels = c("Gulf of Lion","Murcia","Faro","Gulf of Gascogne"),
                    values = brewer.pal(n = 4, name = "RdBu"))+
  ylab("Snps Het/Hom")+
  xlab("Species")+
  theme_classic()+
  scale_colour_manual(name = "Location",
                      labels = c("Gulf of Lion","Costa Calida","Algarve","Bay of Biscay"),
                      values = brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position="bottom")

ppp<-ggplotly(ppp,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
```

### Depth distribution
```{r}
VCF_DEPTH_QUALITY_PER_SITE = VCF_DEPTH_QUALITY_PER_SITE[-1,]
for (i in 1:nrow(VCF_DEPTH_QUALITY_PER_SITE)){
  if (VCF_DEPTH_QUALITY_PER_SITE[i,]$FILTER != "PASS"){
    VCF_DEPTH_QUALITY_PER_SITE[i,]$FILTER = "FILTERED"
  }
}
```

```{r, fig.height=150}
ppp<-ggplot(VCF_DEPTH_QUALITY_PER_SITE[VCF_DEPTH_QUALITY_PER_SITE$DEPTH<=100,], aes(x=DEPTH,fill=FILTER)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.5,binwidth = 1)+
  geom_density(alpha=0.6)+
  geom_vline(xintercept=20,
             linetype="dashed")+
  xlim(c(0,65))+
  scale_fill_viridis_d(begin=0.1,end=0.9)+
  labs(title="Depth",x="Depth", y = "Freq.")+
  theme_classic()+
  facet_grid(SPECIES ~ .)

ppp<-ggplotly(ppp,height=2500)  %>% partial_bundle() %>% hide_legend()
ppp
```

### Quality distribution

```{r, fig.height=150}
ppp<-ggplot(VCF_DEPTH_QUALITY_PER_SITE, aes(x=QUALITY,fill=FILTER)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.5,binwidth = 100)+
  geom_density(alpha=0.6)+
  geom_vline(xintercept=20,
             linetype="dashed")+
  xlim(c(0,5e3))+
  scale_fill_viridis_d(begin=0.1,end=0.9)+
  labs(title="Quality",x="log[Quality]", y = "Freq.")+
  theme_classic()+
  facet_grid(SPECIES ~ .)

ppp<-ggplotly(ppp,height=2500)  %>% partial_bundle() %>% hide_legend()
ppp
```

Packages used {data-icon="fa-map"}
=======================================================================

```{r}
installed.packages()[names(sessionInfo()$otherPkgs), "Version"]
```